/*
 *   This program is free software: you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation, either version 3 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 *    RandomVariates.java
 *    Copyright (C) 2002-2012 University of Waikato, Hamilton, New Zealand
 *
 */

package weka.core;

import java.util.Random;

/**
 * Class implementing some simple random variates generator.
 *
 * @author Xin Xu (xx5@cs.waikato.ac.nz)
 * @version $Revision$
 */
public final class RandomVariates extends Random {

    /** for serialization */
    private static final long serialVersionUID = -4763742718209460354L;

    /**
     * Simply the constructor of super class
     */
    public RandomVariates() {
        super();
    }

    /**
     * Simply the constructor of super class
     *
     * @param seed the seed in this random object
     */
    public RandomVariates(long seed) {
        super(seed);
    }

    /**
     * Simply use the method of the super class
     *
     * @param bits - random bits
     * @return the next pseudorandom value from this random number generator's
     *         sequence.
     */
    protected int next(int bits) {
        return super.next(bits);
    }

    /**
     * Generate a value of a variate following standard exponential distribution
     * using simple inverse method.
     * <p>
     *
     * Variates related to standard Exponential can be generated using simple
     * transformations.
     * 
     * @return a value of the variate
     */
    public double nextExponential() {
        return -Math.log(1.0 - super.nextDouble());
    }

    /**
     * Generate a value of a variate following standard Erlang distribution. It can
     * be used when the shape parameter is an integer and not too large, say, <100.
     * When the parameter is not an integer (which is often named Gamma
     * distribution) or is large, use <code>nextGamma(double a)</code> instead.
     *
     * @param a the shape parameter, must be no less than 1
     * @return a value of the variate
     * @exception Exception if parameter less than 1
     */
    public double nextErlang(int a) throws Exception {
        if (a < 1)
            throw new Exception("Shape parameter of Erlang distribution must be greater than 1!");

        double product = 1.0;
        for (int i = 1; i <= a; i++)
            product *= super.nextDouble();

        return -Math.log(product);
    }

    /**
     * Generate a value of a variate following standard Gamma distribution with
     * shape parameter a.
     * <p>
     * If a>1, it uses a rejection method developed by Minh(1988)"Generating Gamma
     * Variates", ACM Trans. on Math. Software, Vol.14, No.3, pp261-266. <br>
     * If a<1, it uses the algorithm "GS" developed by Ahrens and Dieter(1974)
     * "COMPUTER METHODS FOR SAMPLING FROM GAMMA, BETA, POISSON AND BINOMIAL
     * DISTRIBUTIONS", COMPUTING, 12 (1974), pp223-246, and further implemented in
     * Fortran by Ahrens, Kohrt and Dieter(1983) "Algorithm 599: sampling from Gamma
     * and Poisson distributions", ACM Trans. on Math. Software, Vol.9 No.2,
     * pp255-257.
     * <p>
     * 
     * Variates related to standard Gamma can be generated using simple
     * transformations.
     *
     * @param a the shape parameter, must be greater than 1
     * @return a value of the variate
     * @exception Exception if parameter not greater than 1
     */
    public double nextGamma(double a) throws Exception {
        if (a <= 0.0)
            throw new Exception("Shape parameter of Gamma distribution" + "must be greater than 0!");
        else if (a == 1.0)
            return nextExponential();
        else if (a < 1.0) {
            double b = 1.0 + Math.exp(-1.0) * a, p, x, condition;
            do {
                p = b * super.nextDouble();
                if (p < 1.0) {
                    x = Math.exp(Math.log(p) / a);
                    condition = x;
                } else {
                    x = -Math.log((b - p) / a);
                    condition = (1.0 - a) * Math.log(x);
                }
            } while (nextExponential() < condition);
            return x;
        } else { // a>1
            double b = a - 1.0, D = Math.sqrt(b), D1, x1, x2, xl, f1, f2, x4, x5, xr, f4, f5, p1, p2, p3, p4;

            // Initialization
            if (a <= 2.0) {
                D1 = b / 2.0;
                x1 = 0.0;
                x2 = D1;
                xl = -1.0;
                f1 = 0.0;
            } else {
                D1 = D - 0.5;
                x2 = b - D1;
                x1 = x2 - D1;
                xl = 1.0 - b / x1;
                f1 = Math.exp(b * Math.log(x1 / b) + 2.0 * D1);
            }

            f2 = Math.exp(b * Math.log(x2 / b) + D1);
            x4 = b + D;
            x5 = x4 + D;
            xr = 1.0 - b / x5;
            f4 = Math.exp(b * Math.log(x4 / b) - D);
            f5 = Math.exp(b * Math.log(x5 / b) - 2.0 * D);
            p1 = 2.0 * f4 * D;
            p2 = 2.0 * f2 * D1 + p1;
            p3 = f5 / xr + p2;
            p4 = -f1 / xl + p3;

            // Generation
            double u, w = Double.MAX_VALUE, x = b, v, xp;
            while (Math.log(w) > (b * Math.log(x / b) + b - x) || x < 0.0) {
                u = super.nextDouble() * p4;
                if (u <= p1) { // step 5-6
                    w = u / D - f4;
                    if (w <= 0.0)
                        return (b + u / f4);
                    if (w <= f5)
                        return (x4 + (w * D) / f5);

                    v = super.nextDouble();
                    x = x4 + v * D;
                    xp = 2.0 * x4 - x;

                    if (w >= f4 + (f4 - 1) * (x - x4) / (x4 - b))
                        return xp;
                    if (w <= f4 + (b / x4 - 1) * f4 * (x - x4))
                        return x;
                    if ((w < 2.0 * f4 - 1.0) || (w < 2.0 * f4 - Math.exp(b * Math.log(xp / b) + b - xp)))
                        continue;
                    return xp;
                } else if (u <= p2) { // step 7-8
                    w = (u - p1) / D1 - f2;
                    if (w <= 0.0)
                        return (b - (u - p1) / f2);
                    if (w <= f1)
                        return (x1 + w * D1 / f1);

                    v = super.nextDouble();
                    x = x1 + v * D1;
                    xp = 2.0 * x2 - x;

                    if (w >= f2 + (f2 - 1) * (x - x2) / (x2 - b))
                        return xp;
                    if (w <= f2 * (x - x1) / D1)
                        return x;
                    if ((w < 2.0 * f2 - 1.0) || (w < 2.0 * f2 - Math.exp(b * Math.log(xp / b) + b - xp)))
                        continue;
                    return xp;
                } else if (u < p3) { // step 9
                    w = super.nextDouble();
                    u = (p3 - u) / (p3 - p2);
                    x = x5 - Math.log(u) / xr;
                    if (w <= (xr * (x5 - x) + 1.0) / u)
                        return x;
                    w = w * f5 * u;
                } else { // step 10
                    w = super.nextDouble();
                    u = (p4 - u) / (p4 - p3);
                    x = x1 - Math.log(u) / xl;
                    if (x < 0.0)
                        continue;
                    if (w <= (xl * (x1 - x) + 1.0) / u)
                        return x;
                    w = w * f1 * u;
                }
            }

            return x;
        }
    }

    /**
     * Main method for testing this class.
     *
     * @param ops # of variates/seed, default is 10/45
     */
    public static void main(String[] ops) {

        int n = Integer.parseInt(ops[0]);
        if (n <= 0)
            n = 10;
        long seed = Long.parseLong(ops[1]);
        if (seed <= 0)
            seed = 45;
        RandomVariates var = new RandomVariates(seed);
        double varb[] = new double[n];

        try {
            System.out.println("Generate " + n + " values with std. exp dist:");
            for (int i = 0; i < n; i++) {
                varb[i] = var.nextExponential();
                System.out.print("[" + i + "] " + varb[i] + ", ");
            }

            System.out.println("\nMean is " + Utils.mean(varb) + ", Variance is " + Utils.variance(varb) + "\n\nGenerate " + n + " values with" + " std. Erlang-5 dist:");

            for (int i = 0; i < n; i++) {
                varb[i] = var.nextErlang(5);
                System.out.print("[" + i + "] " + varb[i] + ", ");
            }

            System.out.println("\nMean is " + Utils.mean(varb) + ", Variance is " + Utils.variance(varb) + "\n\nGenerate " + n + " values with" + " std. Gamma(4.5) dist:");

            for (int i = 0; i < n; i++) {
                varb[i] = var.nextGamma(4.5);
                System.out.print("[" + i + "] " + varb[i] + ", ");
            }

            System.out.println("\nMean is " + Utils.mean(varb) + ", Variance is " + Utils.variance(varb) + "\n\nGenerate " + n + " values with" + " std. Gamma(0.5) dist:");

            for (int i = 0; i < n; i++) {
                varb[i] = var.nextGamma(0.5);
                System.out.print("[" + i + "] " + varb[i] + ", ");
            }

            System.out.println("\nMean is " + Utils.mean(varb) + ", Variance is " + Utils.variance(varb) + "\n\nGenerate " + n + " values with" + " std. Gaussian(5, 2) dist:");

            for (int i = 0; i < n; i++) {
                varb[i] = var.nextGaussian() * 2.0 + 5.0;
                System.out.print("[" + i + "] " + varb[i] + ", ");
            }
            System.out.println("\nMean is " + Utils.mean(varb) + ", Variance is " + Utils.variance(varb) + "\n");

        } catch (Exception e) {
            e.printStackTrace();
        }
    }
}
